Diffusion entropy and waiting time statistics of hard x-ray solar flares 
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We analyze the waiting time distribution of time distances r between two nearest-neighbor flares. 
This analysis is based on the joint use of two distinct techniques. The first is the direct evaluation 
of the distribution function tp(r), or of the probability, ^(tau), that no time distance smaller than 
a given r is found. We adopt the paradigm of the inverse power law behavior, and we focus on 
the determination of the inverse power index /i, without ruling out different asymptotic properties 
that might be revealed, at larger scales, with the help of richer statistics. The second technique, 
called Diffusion Entropy (DE) method, rests on the evaluation of the entropy of the diffusion process 
generated by the time series. The details of the diffusion process depend on three different walking 
rules, which determine the form and the time duration of the transition to the scaling regime, as 
well as the scaling parameter 5. With the first two rules the information contained in the time 
series is transmitted, to a great extent, to the transition, as well as to the scaling regime. The same 
information is essentially conveyed, by using the third rules, into the scaling regime, which, in fact, 
emerges very quickly after a fast transition process. We show that the significant information hidden 
within the time series concerns memory induced by the solar cycle, as well as the power index \x. 
The scaling parameter S becomes a simple function of fi, when memory is annihilated. Thus, the 
three walking rules yield a unique and precise value of /i if the memory is wisely taken under control, 
or cancelled by shuffling the data. All this makes compelling the conclusion that \i = 2.138 ± 0.01. 
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I. INTRODUCTION 

The study of solar flares is becoming popular among 
the researchers working at the frontier of statistical me- 
chanics, due to the widely shared conviction that they 
are a signature of a significant departure from the con- 
dition of ordinary Brownian motion Gl-0]. As pointed 
out by Wheatland, [|| , the distribution of times between 
flares, gives information on how to model flare statis- 
tics. In this paper we shall be referring to these times, 
denoted by us with the symbol r, as the time distance 
between two consecutive events and the corresponding 
distribution density will be denoted by ip(r). Although 
the agreement on the fact that flare statistics depart from 
ordinary statistical mechanics is general, there seems to 
be the still unsettled issue of what is the proper model 
which will account for this form of anomalous statistics. 
Does this form of statistics reflect self-organized critical- 
ity or turbulence ||? We think that the settlement of 
this delicate issue is made difficult by the fact that, al- 
though many authors claim that ip( T ) is an inverse power 
law with power index \x, the actual value of /i still seems 
to be uncertain. In fact, the authors of Ref. |Ej pro- 
pose fj, = 1.7 and those of Ref. claim that fi = 2 is 
the proper power law index. Boffctta et al. || propose 
/i = 2.4. Finally, Wheatland explains the origin of the 
power law behavior with a model yielding fi = 3.0, [0. 

As it will be made clear by the theoretical analysis of 
this paper, it is possible to prove, without taking posi- 
tion on the origin of the inverse power law behavior, that 
fi = 3 and /i = 2 are critical values. In fact, we note that 



moving from [i > 3 to /i < 3 is equivalent to a phase tran- 
sition from the Gaussian to the Levy basin of attraction 
Q , and moving from /i > 2 to fi < 2 implies a transition 
from the condition of Levy statistics to a form of out of 
equilibrium regime 0. Thus, an uncertainty larger than 
the distance of the border fi = 2 from the border /.t = 3 
is judged by us to be an unsatisfactory condition that 
might delay the settlement of the issues concerning the 
complex dynamics underlying the waiting time statistics. 
The main purpose of this paper is to illustrate a statisti- 
cal method of analysis that yields a reliable value for the 
power index \x. We hope that this result might be useful 
for the researchers in this interesting field of investigation 
and at the same time might be beneficial in general for 
all those who are interested in the statistical analysis of 
time series. 

The outline of the paper is as follows. In Section II 
we review the method of Diffusion Entropy (DE) that 
will be a crucial step of the statistical analysis done in 
this paper. Although the method has been applied some- 
where else [pj-p^l , we will present a short review in order 
to make this paper as self-contained as possible. In Sec- 
tion III, we illustrate a dynamical model that in general 
results in time sequences that are statistically equivalent 
to those observed in real data. This model is not limited 
to the case of inverse power laws, but here we make the 
assumption that the shifted inverse power law is an ideal 
condition convenient to analyze solar flares, and we study 
the explicit form emerging from this condition. In Sec- 
tion IV, we illustrate two walking prescriptions that will 
be used to convert the real data into random trajectories. 
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The benefit of adopting several walking prescriptions was 
discussed in Ref. |£j] . Here we introduce two new rules and 
we apply both of them as well as one of those introduced 
in Ref. ||. In Section V we prove that the numerical 
evaluation of the probability of getting a waiting time 
larger than a given r yields a value for \x more accurate 
than that afforded by the waiting-time distribution iP(t). 
In Section VI we show how to process the data to make 
an efficient use of the DE method. In Section VII we use 
the DE method to further reduce the error of Section V. 
We devote Section VIII to concluding remarks. 



II. DIFFUSION ENTROPY 

The main idea of this approach to scaling is remark- 
ably simple. Let us consider a sequence of M numbers, 
, with i = 1, . . . , M. The purpose of the DE algo- 
rithm is to establish the possible existence of a scaling, 
either normal or anomalous, in the most efficient way as 
possible without altering the data with any form of de- 
trending. Let us select first of all an integer number I, 
fitting the condition 1 < I < M. This integer number 
will be referred to by us as "time" . For any given time I 
we can find M — I + 1 sub-sequences defined by 



i+s j 



0. 



,M — I. 



(1) 



For any of these sub-sequences we build up a diffusion 
trajectory, labeled with the index s, defined by the posi- 
tion 



(0 = E^ (s) = Ee 



z+s ■ 



(2) 



Let us imagine this position as referring to a Brow- 
nian particle that at regular intervals of time has been 
jumping forward or backward according to the prescrip- 
tion of the corresponding sub-sequence of Eq.([j]). This 
means that the particle before reaching the position that 
it holds at time I has been making I jumps. The jump 
made at the «-th step has the intensity \Q \ and is for- 

(s) 

ward or backward according to whether the number Q 
is positive or negative. 

We are now ready to evaluate the entropy of this diffu- 
sion process. To do that we have to partition the x-axis 
into cells of size e(l). When this partition is made we 
have to label the cells. We count how many particles are 
found in the same cell at a given time /. We denote this 
number by Ni(l). Then we use this number to determine 
the probability that a particle can be found in the i-th 
cell at time I, Pi(l), by means of 



Mi) 



_ Ni(l) 



(M -1 + 1)' 



(3) 



At this stage the entropy of the diffusion process at time 
I is determined and reads 



£K0 = -I>(0J»[Pi(0]- 



(4) 



The easiest way to proceed with the choice of the cell size, 
e(l), is to assume it to be independent of / and determined 
by a suitable fraction of the square root of the variance 
of the fluctuation £(£). In the case in which the numbers 
£t are +1, and -1, e = 1 is the natural choice. 

Before proceeding with the illustration of how the DE 
method works, it is worth making a comment on how to 
define the trajectories. The method we are adopting is 
based on the idea of a moving window of size / that makes 
the s — th trajectory closely correlated to the next, the 
(s + 1) — th trajectory. The two trajectories have I — 1 
values in common. It is worth making a comparison with 
the technique of Detrended Fluctuation Analysis (DFA) 
PJ . The DFA is a popular method of scaling analysis, 
aiming at detecting the long-range correlations in seem- 
ingly non-stationary time series that in the last few years 
has been used in more than 100 publications The 
DFA is based on non-overlapping windows, and, conse- 
quently, trajectories with different labels are totally in- 
dependent from one another. The motivation for using 
overlapping windows, with the DE method, is given by 
our wish to establish a connection with the Kolmogorov- 
Sinai (KS) entropy jl3],[l4|]. In Section III we shall make 
further comments on this connection. The KS entropy 
of a symbolic sequence is evaluated by moving a window 
of size I along the sequence. Any window position cor- 
responds to a given combination of symbols, and from 
the frequency of each combination it is possible to de- 
rive the Shannon entropy S(l). The KS entropy is given 
by the asymptotic limit limi—tooSQ,)/!. We believe that 
the same sequence, analyzed with the DE method, at the 
large values of I where a finite KS entropy shows up, must 
yield a well defined scaling 5. To realize this correspon- 
dence we carry out the determination of the DE by using 
the same criterion of overlapping windows as that behind 
the KS entropy. 

Details on how to deal with the transition from the 
short-time regime, sensitive to the discrete nature of the 
process under study, to the long-time limit where both 
space and time can be perceived as continuous, are given 
in Ref. (To) . Here we make the simplifying assumption of 
considering large enough times as to make the continuous 
assumption valid. In this case, the trajectories, built up 
with the above illustrated procedure, correspond to the 
following equation of motion: 



dx 
It 



(•5) 



where denotes the value that the time series under 
study gets at the I — th site of the sequence under study. 
This means that the time t = I (with I 3> 1) is thought of 
as a continuous and that the function £(Z) is a function of 
this continuous time. In this case the Shannon entropy 
reads 
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S(t) = - dxp(x,t)ln\p(x,t)]. (6) 

J — OO 

We also assume that 

p^ t ) = j^m F (jk) (7) 

and that F(y) maintains its form, namely that the statis- 
tics of the process is independent of time. Let us plug 
Eq.(^) into Eq. (||). Using a simple algebra, we get: 



where 



and 



S(t) = A + S(t)t, 

I dyF(y)ln[F(y)} 

) 

= ln(t). 



(8) 
(9) 
(10) 



The assumptions made to get the result of Eq. (j|) 
are not correct during the transition process, and conse- 
quently the DE method can be used as a reliable way to 
detect scaling only in the long-time limit. The DE can 
be used however to shed light into the regime of transi- 
tion that is deeply connected with the foundation itself 
of statistical mechanics. According to Khinchin Jl5| the 
central limit theorem is fundamental for the realization 
of canonical equilibrium. As well known, a process re- 
sulting from the sum of N independent variables yields 
a Gaussian distribution, provided that N is large and 
the single variables have a probability distribution with 
a finite second moment. A physical process making N 
increase from values of the order of unity to values so 
large as to fit the prediction of the central limit theorem 
can be perceived as a transition from the microscopic to 
the macroscopic regime, where thermodynamics applies. 
If the microscopic variables do not have a finite second 
moment, the ordinary central limit theorem must be re- 
placed by the generalized central limit theorem [|l6| and 
in the limiting case of N — > oo we find Levy rather than 
Gauss statistics. We can generalize the point of view of 
Khinchin and consider also in this case the process of 
transition of N from small to large values as a form of 
transition from the microscopic to the thermodynamic 
regime. 

Due to the nature of the DE method, the role of N 
is here played by the "time" t. The microscopic regime 
refers to the fluctuation of £j and the macroscopic regime 
corresponds to the fluctuations of the diffusion coordi- 
nate x{t). The time evolution of S(t) towards the final 
value, independent of time, reflects the transition from 
dynamics to thermodynamics. 

We shall adopt three different walking rules (see Sec- 
tion IV). The first two rules are characterized by an ex- 
tended regime of transition from dynamics to thermo- 
dynamics. Notice that the real data available are fi- 
nite, thereby producing saturation effects in the long- 
time regime. Consequently, the region where the ideal 



scaling shows up, is an intermediate time region follow- 
ing the extended initial transition and preceding the long- 
time saturation regime. This has the effect of reducing 
the size of the time region that can be fruitfully used for 
scaling detection. As we shall see, this is the reason why 
the DE method must be supplemented by the use of ar- 
tificial sequences. The third rule, on the contrary, yields 
a fast transition to the thermodynamic regime, thereby 
allowing us to determine the scaling by a direct use of 
Eq. (|J), with S(t) assuming the time independent value 
of the thermodynamic limit. 



III. DYNAMIC MODEL 

The solar flares analyzed in this paper are perceived as 
a sequence of events occurring at unpredictable times, U, 
with i = 1, . . . , M, where M is the label of the last event 
considered. We do not take into account the intensity 
of these events, which will be studied somewhere else. 
Thus, the most important property for us to study, is the 
time distribution density, ip{r), with r denoting the time 
distance between two nearest neighbor events, Tj = — 
tj. Let us make the assumption that the experimental 
analysis of the time series yields the form 



1>(t) = (m - 1) 



rpfl- 



(T + T)f 



(11) 



We make the key assumption that the numbers t% are 
uncorrelated. As we shall see, the theory of this paper 
affords also a criterion to assess if this crucial assumption 
is correct or not. Under this key assumption we can build 
up a dynamic model that is statistically equivalent to the 
solar dynamics generating the sequence of the U's. Let 
us consider the dynamic process: 



dy/dt = \y z 



(12) 



with z > 1. Let us imagine that the trajectory y(t) moves 
within the interval [0, 1]. Let us assume also that when 
the trajectory reaches the right border of this interval it 
is injected back within this interval by means of a random 
selection of the initial position y(0). The random selec- 
tion is done by using a random number generator that 
assigns the same probability to the numbers of the inter- 
val [0,1]. The connection between the initial condition 
and the exit time r is given by 



y(0) = [! + (*- l)Ar] 



(13) 



This leads immediately to the distribution of Eq.(|l] 
with 



z 



and 



T = 



z-1 



M-l 



(14) 



(15) 



3 



We are now equipped to establish a connection with the 
entropy production per unit of time. Randomness here is 
involved at the moment of selecting the initial condition, 
and is characterized by an unknown amount of entropy 
increase, H. If /j, > 2, the distribution of Eq.(|Tl|) yields 
a finite mean waiting time 



< T > = 



T 



(16) 



It is evident then that the rate of entropy production per 
unit of time is given by 



h E = H 



T 



(17) 



Using the dynamic model it is possible to establish a more 
proper connection with the KS entropy. However, this is 
of no great relevance within the context of the present 
paper. Therefore, we limit ourselves to considering the 
entropy production of Eq. ([l7]), where the subscript E 
stands for "external". In fact, in the picture adopted 
in the present paper the source of entropy production is 
the random selection of the numbers of the interval [0,1], 
an action external to the process under study. It has to 
be pointed out that this external entropy production is 
subtly related to the KS entropy, which, on the contrary, 
is interpreted as being of internal origin p7y . This is 
so because the dynamical model is a map with a very 
sharp chaotic region that reduces to a set of zero measure, 
confined to the point y — 1, in the limiting condition 
where the idealized model of this section applies. 

In the case fi < 2 the entropy produced is proven |?]] 
to be the following function of time: 



S{t) oc f" 



(18) 



It is evident that in the limiting case of very large time 
values the entropy production per unit of time vanishes, 
thereby implying that the condition, /i = 2 is a border 
at which a kind of phase transition occurs. In the region 
fi > 2 the dynamical system of Eq. ( |l2| ) has an invari- 
ant distribution. In the region /x < 2 the system does 
not have an invariant distribution [Q. From an intuitive 
point of view we can imagine that during the observation 
process the system keeps moving towards an equilibrium 
distribution, as a kind of Dirac delta function located at 
V = @- The time necessary to reach this invariant 
distribution is infinite. 

In conclusion, an infinitesimally small change from 
ji > 2 to ji < 2 would have the effect of annihilating the 
invariant distribution and of making the process "non- 
stationary." The method of analysis of this paper will 
allow us to assess that = 2.138 ±0.01, namely, that the 
solar flares fluctuations are stationary, even if very close 
to the border with the "non-stationary" region. This 
result will be obtained by a direct evaluation of /i, sup- 
plemented by the adoption of the DE method. As we 
shall see, this conclusion is reached after settling a ma- 
jor problem caused by the existence of a genuine form 



of non-stationary behavior, a fact that is more properly 
related to dynamic rules changing upon change of time. 
This will lead us to the final conclusion that the model 
of Eq. (12) is a fairly accurate way of mimicking solar 



flare dynamics, with z < 2 (fi > 2). In Sections VII and 
VIII we shall make some conjectures on how to improve 
this model to take into account the time dependence of 
the solar flare rate. 



IV. ON THREE DISTINCT PRESCRIPTIONS TO 
WALK 

The scaling detected by the DE method is not indepen- 
dent of the walking rules that we adopt. The outcomes 
of DE method are not unique, due to the dependence of 
the scaling parameter S on the walking rules, and this 
casts doubts on this method of analysis. However, the 
task of this analysis is an indirect evaluation of the wait- 
ing time distribution ip(r), or, equivalently, in the inverse 
power law case, of the index fi. If we take for granted the 
inverse power law structure of ip( T )> the- power index /x 
is unique. We adopt the following prescriptions for the 
random walker: 

Rule No. 1. Make a jump of fixed intensity, only when 
you meet an event, and do it always in the same direction. 

Rule No. 2. As with Rule No. 1, make a jump only 
when you meet an event, but do it either in the positive 
or negative direction according to a coin tossing prescrip- 
tion. 

Rule No. 3. Walk at fixed interval of times, with 
jumps in the same direction, of intensity proportional to 
the time distance between two nearest-neighbor events. 

Note that here we analyze the sequence {t^}, where 
each value denotes the time distance between two 
nearest-neighbor flares (regarded as events). Thus, rules 
No. 1 and No. 2 imply that the random walker makes in- 
stantaneous jumps at the times of flare occurrence. With 
rule No. 3 the random walker, at times t = 1, 2, n, 
makes jumps ahead of intensity equal to the values Ti of 
the sequence under study. Note that rule No. 1 is one of 
the two rules used in Ref. fi. Here we use for the first 
time rules No. 2 and No. 3. Using the theory of Ref. 

which, in turn, essentially rests on the generalized 
central limit theorem |l6| and on the work of Feller [|l8| , 
we obtain the following prescriptions: 




1). 



5 = 



0.5, 



1 < n < 2 

2 < n < 3 
H > 3, 



1< /i < 2 
/i > 2 



and 



S = l/(ji-l), M>U 



(19) 



(20) 



(21) 
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for rules No. 1, No. 2 and No. 3, respectively. 

Fig. 1 shows clearly that the adoption of rule No. 1 
alone would yield two distinct possible values for \i when 
S gets values within the interval [0.5,1]. However, the 
joint adoption of this and the other two rules settles this 
ambiguity. We also notice that both rule No. 1 and 
rule No. 2 reflect the phase-transition character of the 
condition fi — 2, while rule No. 3, apparently, does not. 
However, we see that rule No. 3 for fx < 2 yields a value of 
S > 1, namely a diffusion process faster than the ballistic 
diffusion. This is a consequence of the non-stationary 
nature of the condition fi < 2. 

It is important to stress that these rules imply that 
the numbers Tj are not correlated. Furthermore, these 
rules rest on the assumption that the asymptotic limit of 
V>( T ) i s an inverse power law distribution with no trun- 
cation. We shall see that the DE method is sensitive to 
the correlation among the numbers Tj, and that the 11- 
year solar cycle is responsible for that correlation. As 
to the truncation of the inverse power law at the large 
distances, this is another delicate issue worth of some 
comments. Laherrere and Sornette, |l9|, suggest that 
the stretched exponential family might have a theoreti- 
cal motivation stronger than the power-law distribution. 
On the other hand, in the intermediate time region a 
stretched exponential is indistinguishable from a power 
law. The two proposed fitting functions become distin- 
guishable one from the other in the long-time regime, 
which is affected by poor statistics. However, the work 
of Refs. and [glj show that a truncation of the power 
law of il>( T ) & t large times yield an ultra-slow conver- 
gence to normal diffusion, with effects that are beyond 
the range of observation of the DE analysis, due to the 
data statistical limitation. 

We shall see that both rule No. 1 and rule No. 2 yield 
a very slow transition to the scaling regime. Due to the 
statistical limitation of our data, the scaling regime turns 
out to be a relatively short time region between transition 
and saturation regime. Thus, we shall be forced to carry 
out our analysis with the help of artificial sequences with 
the same number of terms as the real data, by fitting the 
DE curves produced by the real data with the DE curves 
generated by the artificial sequences. The adoption of 
the third rule, on the contrary, yield a fast transition 
to the thermodynamic regime and, consequently, makes 
it possible to make a direct evaluation of 6. In both 
cases, however, the physical consequences of a possible 
truncation of the inverse power law are beyond our range 
of observation. 



V. STATISTICAL ANALYSIS OF THE REAL 
DATA: V(t) AND *(r) 

In this section we plan to derive the waiting time dis- 
tribution ip(r) directly from the statistical analysis of the 
real data, the x-rays emitted by solar flares in the case 



here under study. At first sight, one might think that a 
direct determination of i/)(t) is more convenient than any 
indirect approach. Actually, it is not so. As mentioned in 
Section I, we find that the evaluation of the probability 
of finding no time distance larger than a given r, denoted 
by ^(t), defined by 

oo 

¥(r) = J mdt, (22) 

r 

is more convenient than the direct evaluation of ip( T )- 
In later sections we shall prove a striking property: the 
evaluation of /i through the DE method, an approach less 
direct than the evaluation of 'J'(t), is still more efficient. 

The data are a set of 7212 hard x-ray peak flaring 
event times obtained from the BATSE/CGRO (Burst 
and Transient Source Experiment aboard the Compton 
Gamma Ray observatory satellite) solar flare catalog list. 
The data is a nine-year series of events from 1991 to 2000. 
If the time At between two consecutive solar flares is ex- 
pressed in seconds, the range goes from 45 to 10,000,000 
seconds, as shown in Fig. 2. Fig. 3 shows the rate of 
solar flares per month from April 1991 to May 2000. The 
set of data studied here concerns a time period of 9 years, 
and, consequently, a large part of the whole 11-year so- 
lar cycle. Fig. 3 shows that during a large portion of 
this 11-year cycle the flare rate undergoes big changes, 
thereby significantly departing from the uniform distribu- 
tion. Furthermore, it is worth remarking that, as shown 
by Fig. 4, the 11-year solar cycle is not a mere harmonic 
oscillation with the period of 11 years, but a complex 
dynamic process with many components. 

The direct evaluation of the waiting time distribution, 
"0(t), needs the data to be distributed over many bins 
with the same size. When only a few data are available, 
the bin size cannot be too small, and, in turn, the adop- 
tion of bins of large size can produce incorrect power law 
indices. In proceeding with the direct evaluation of the 
key parameter /z, first of all, we have to adopt a proper 
criterion to determine the size A, of the i-th bin. We note 
that the waiting time distribution is expected to be an 
inverse power law. If we adopted bins of equal size, those 
corresponding to large times would collect a very lim- 
ited amount of data, thereby resulting in a non reliable 
evaluation of the frequencies. To bypass this difficulty we 
adopt bin sizes that are constant in the logarithmic scale. 
This means that ln{ji) — £n(Tj_i), where Tj and Tj_i are 
the middle times of two consecutive bins, is constant. We 
define the width of the i-th bin as A = r, — t,_i, thereby 
making it become an exponentially increasing function of 
the sequence position, so as to widely compensate for the 
density decrease. In this representation the probability 
density ipfc) is expressed by 

= ^ , (23) 
where N is the total number of data points, Nj is number 
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of points located within the i-th bin, and Aj, as earlier 
said, is the width of the i-th bin. 

The fitting is done by using the prescription of a power 
law of the type of Eq. fll] 



(T- 



(24) 



with A\ , T and fj, being three independent fitting param- 
eters. It is worth noting that the normalization condi- 
tion reduces the three independent parameters to two, as 
made clear by Eq. ( pd| ) , which is a function of only T and 
/x. We find it to be more convenient to adopt three inde- 
pendent fitting parameters, with the understood proviso 
that the departure of A\ from the value (/i — 1)T' I_1 can 
be interpreted as a way to estimate the inaccuracy of the 
adopted fitting procedure. 

The fitting is done by using an implementation of the 
nonlinear least-squares (NLLS) Marquardt-Levenberg al- 
gorithm p2[ . The NLLS algorithm may not give unique 
values for the fitting parameters. It needs initial guesses 
for the free parameters and the final results may change 
or be affected by huge errors. This fitting procedure 
yields: T = 8787, \i = 2.12 ± 0.32 and A x = 31006. The 
evaluated value of A\ is not far from the value 29236 that 
would be required by the normalization condition. How- 
ever, there are very large errors of the order of 100%, with 
an error on the parameter fj, of the order of 15%, thereby 
implying 1.80 < /i < 2.44. This means that the result 
of this fitting procedure would prevent us from assessing 
the important question raised in Section III on whether 
the process is stationary or non stationary. The large 
error of this procedure depends upon the initial values 
assigned to the three fitting parameters, T, [i and A\, 
whose choice requires a more efficient criterion. It also 
depends on the fact that there are oscillations around the 
fitting curve, as clearly illustrated by Fig. 5. 

As earlier mentioned several times, a more accurate 
fitting is obtained using the function ^(r). Again we do 
not pay attention to the normalization constraints and 
we adopt the following fitting function 



#(t) = A 2 



1 



T 



(25) 



As shown by Fig. 6, the fitting of the real data is now 
much more accurate than that of Fig. 5. The fitting pa- 
rameters used are: A 2 = 30657 ± 16590, T = 8422 ±500, 
fi = 2.144 ± 0.05. This sets on the key parameter /i the 
constraint 2.094 < // < 2.194, which has the very attrac- 
tive property of establishing the stationary nature of the 
dynamic model behind the solar flares fluctuations. The 
results of this search for fi, based on the direct evaluation 
of iP(t) and on the use of ^(t), are summarized in Table 
I. We note that the uncertainty interval associated with 
the use of 'J'(t) is contained within the wider uncertainty 
interval produced by the use of i/j(t). This means that 
we are coming closer to the real value of [i. The width of 



the uncertainty interval will be further reduced by using 
the DE method. 

Table I 







1.80 < \l < 2.44 


2.094 <fi< 2.194 



VI. DIFFUSION ENTROPY OF SOLAR FLARES. 

This section is devoted to the analysis of the solar flares 
data by means of the DE method. The final result will 
be given by /i = 2.138 ± 0.01, namely a value for /i even 
more accurate than that obtained in Section V by using 
^(t). Wc shall prove also that the DE method allows us 
to establish some aspects of the dynamics behind solar 
flares that would be overlooked by an analysis based only 
on the use of the waiting time distribution. 

The first issue that we have to solve is how to process 
the data so as to apply the three walking rules of Section 
IV. The data accessible to us are the times Tj = t{ — ij-i, 
with ti and tj_i denoting the time of occurrence of the i- 
th and the (i-l)-th solar flare, respectively. However, the 
direct adoption of these numbers would result in technical 
difficulties that are bypassed by referring ourselves to the 
new sequence of numbers 



[3j = Int 



At* 
A 



1. 



(26) 



where i~ni[x] denotes the integer part of x. The adop- 
tion of A = 1 would be virtually equivalent to referring 
ourselves to the original sequence of numbers. However, 
preliminary trials with changing values of A led us to 
conclude that there are problems with the adoption of 
both excessively small and excessively large values of A. 
The adoption of excessively small values of A would make 
the computer analysis too slow and would require an ex- 
cessively large amount of computer memory. This is the 
reason why we cannot use the original sequence of num- 
bers. The adoption of excessively large values of A, on 
the other hand, would produce statistical saturation, and 
a consequent sub-diffusion process that would not accu- 
rately reflect the dynamics behind the data. We adopted 
the criterion of using the largest value of A compatible 
with negligible saturation effect. Preliminary attempts 
made it possible for us to assess that this convenient value 
is given by A = 3600. 

After processing the data, we have to realize the three 
walking rules of Section IV. We note that according to 
the prescription of Section II, diffusion is generated by 
the random walker jumping at any time step. The ran- 
dom walker makes jumps of intensity ahead or back- 
ward, according to whether ^ > or ^ < 0. Thus, we 
create a new sequence of 0's and l's, with the fol- 
lowing prescription. We consider a sequence of infinite 
empty sites, labeled by the integer index i, considered as 
a discrete time, running from i = 1 to i = oo. We divide 
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this sequence into patches of width j3j . The first patch 
consists of the sites i — 1, i — 2,.., i — f3i, the second 
patch consists of the sites i = 0i + 1, (3i + 2, .., 0i+ fe, 
and so on. We assign the value to all the sites of the 
same patch but the last site. This means that the random 
walker walks only at the end of the patch, namely, at the 
occurrence time of an event. To apply rule No. 1, with 
the random walker always moving in the same direction, 
we always assign to the last site of a given patch the value 
of 1. To apply rule No. 2 we assign to the last site of any 
patch either the value 1 or the value -1, according to the 
coin tossing rule. The coin tossing prescription is realized 
by using a random number generator. To reduce the risk 
of artificial periodicity we create 10 different sequences, 
each corresponding to a different random distribution of 
l's and -l's. For any sequence we run the DE method 
and then we make the average over the 10 resulting DE 
curves. To apply the rule No. 2, which will be shown in 
action in Section VII C, we have to identify with /3,-. 

The DE results obtained applying rule No. 1 are illus- 
trated in Fig. 7. This figure shows one of the benefits of 
the DE method. According to rule No. 1, we have to use 
the prescription of Eq. fll9| ) . The most accurate of the 
values of fx, discussed in Section V, is fx = 2.144. This 
value, being smaller than 3 and larger than 2, makes us 
adopt the formula 5 = 1/(11— 1), and yields the scaling 
parameter S — 0.874, which is the slope of the straight 
line of Fig. 7. 

This theoretical prediction implies that the times Tj of 
the sequence {7;} are not correlated with each other. In 
the specific case of seasonal periodicity described by har- 
monic oscillations, the numerical results of Ref. || prove 
that the scaling detected by the DE, as well as by other 
methods to detect scaling, is higher than the Brownian 
motion scaling 8 = 0.5. This is so even when there is no 
correlation in addition to seasonal periodicity. We elim- 
inate this effect, by shuffling the data. The DE method 
can be applied to both the original sequence of 0i and 
to the shuffled sequence. If the DE yields two different 
curves, this is a proof of the fact that there is memory 
in the original sequence. This is an important property 
that cannot be revealed by the analysis of the waiting 
time distribution, t/>(t). Fig. 7 shows that this is the 
case. In fact we see that the DE curve corresponding to 
the shuffled data, after the transition region at short time 
and before saturation, has a slope distinctly smaller than 
the curve referring to the non shuffled data. Further- 
more, this slope is closer to the slope of the straight line 
corresponding to the finding of Section V, which yields 
ll = 2.144, and, consequently, according to Eq. (19), 
S = 0.874. However, both shuffled and non-shuffled data 
yield saturation effects at a time scale of the order of 
tsat = 1, 500 hours. These saturation effects set limits to 
the accuracy of the determination of the value of /1 by 
means of the DE method. 

In Fig. 8 we illustrate the results obtained by using rule 
No. 2. It is remarkable that in this case the shuffled data 
yield, with the DE method, an entropy increase faster 



(rather than slower) than the non-shuffled data. This is 
a consequence of the fact that in this case the deviation 
from ordinary diffusion, produced by time periodicity, 
would generate sub-diffusion rather than super-diffusion. 
We notice that the difference between the shuffled and 
non-shuffled curves is smaller than that in the case of 
Fig. 7 (rule No. 1) and that the saturation effects show 
up at later times. We thus conclude that rule No. 2 
is much less sensitive to periodicities and to saturation 
effects than rule No. 1. 



VII. A FURTHER IMPROVEMENT: USE OF 
ARTIFICIAL SEQUENCES 

We have seen that the DE method reveals the exis- 
tence of memory effects that are overlooked by the direct 
evaluation of the waiting time distribution. However, as 
pointed out in Section II and illustrated by the numeri- 
cal results of Section VI, the time region where the DE 
method might be fruitfully used to detect scaling, is re- 
duced to an intermediate time region, after the transition 
from dynamics to thermodynamics, and before the sat- 
uration effects. This has the unwanted effect of setting 
limitations to the accuracy of the DE method. To by- 
pass this difficulty we generate artificial sequences with 
the same statistical limitations of the real data, and then 
we search for the parameter fi that establishes the most 
accurate fitting with the DE curves derived from real 
data. 

To make this procedure as reliable as possible we pro- 
ceed as follows. We assume that ip(r) has the form 



4>(t) = 



A 



(t + T y 



(27) 



where T and /i are our fitting parameters. The constant 
A is determined by the normalization condition through 



A 



1 



45 (T + rY 



-dr 



(28) 



The fitting parameters are made to change around the 
mean values established by the results of Section V which 
yield fi = 2.144 ± 0.05 and T = 8422 ± 500. Note that 
in the real data no time exists with a value smaller than 
t = 45 sec. This is the reason why the integration in 
Eq. (28) is done from 45 to 00 rather than from to 00. 
The number of data available to us are 7211. Thus we 
produce 7211 values of Tj, according to the prescription 



1 



(T + 45) 



VL-l 



it£ - 1) m 

A 



T 



(29) 



with the number j/j randomly selected in the interval 
[0,1]. It is straightforward to prove that the resulting 
distribution of Tj is the same as that of Eq. ( p7|) and fits 
the condition of Eq. (Eq) . At this stage we are ready to 
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compare the DE curves generated by the artificial data 
to the DE curves generated by the real data, using both 
rule No. 1 and rule No. 2. The comparison is made 
with the DE curves corresponding to shuffled data, since 
the artificial sequences arc generated without correlation 
among the numbers 7$. 

Let us discuss first the results concerning rule No. 1. 
These results are illustrated in Figs. 9. In Fig. 9a we 
show the effect of changing /i in the interval [2.094, 2.194] , 
with T = 8422 and in Fig. 9b we show the effect of chang- 
ing T in the interval [7922, 8922], with \i = 2.144. We see 
that the DE curves of the artificial sequences fluctuate 
within an error strip containing the DE curve of the real 
data. The size of this error strip increases upon change of 
time and we see that the spreading caused by the change 
of T is much smaller than that caused by the change of /x. 
From a qualitative point of view, the results concerning 
rule No. 2, shown in Figs. 10a and 10b, are very similar. 



A. A more accurate measurement of fi. 

We have seen that the area of the T-error strip is sig- 
nificantly smaller than that of the /u-error strip, at least 
five times smaller. Therefore, we can improve the accu- 
racy of fj, by assigning to T a fixed value and looking for 
the value of fi ensuring the best fitting of the real data. 
We assign to T the value of 8422, and we proceed with 
the search for the best fitting. The results are illustrated 
in Figs. 11a and lib. The result concerning rule No. 
1 is good, as seen in Fig. 11a. As expected, Fig. lib 
shows that the result concerning rule No. 2 is even bet- 
ter, and we think that it can be judged to be excellent. 
This extremely accurate result is due to the DE curve of 
the artificial sequence coinciding with the DE curve of 
real data over the wide range of 1000 hours of diffusion. 
On the basis of this excellent fitting, we conclude that 



/! = 2.138 ±0.01 



(30) 



B. Non shuffled data and an artificial sequence with 
suitable memory. 



In Section VI, we have noticed that the result of the 
DE analysis depends on whether the real data are shuf- 
fled or not. We think that in the original data there are 
signs of the 11-year solar cycle and other subcycles. This 
makes it harder to establish a connection between the 
scaling S and the power index ji. However, if our conclu- 
sion that /i = 2. 138 ±0.01 is correct, it should be possible 
to fit the DE curve of the non-shuffled original data with 
no further change of the fitting parameters T and //, pro- 
vided that we sort the artificial sequence in such a way 
as to mimic the solar periodicity. Rather than doing that 
with a model, for instance a suitable modulation of the 
parameter A of Eq. (12), we proceed in a more direct 



way, according to the following procedure. Let us call 
Ri and Ai the i — th numbers of the real and artificial 
sequence used in subsection A, respectively. The i — th 
number of the sorted artificial sequence is denoted by Si. 
The subscript i ranges from 1 to N. The number Si is 
fixed by selecting from the set of A^s the number that 
is closest to Ri, this being, let us say, Aj^y We thus 
set Si = ■ The number Ajm is eliminated from 
the artificial sequence. Then, we move to i?2 and from 
the set of the remaining N-l numbers of the artificial 
sequence we select the closest one to it, this being, let 
us say, Aj( 2 )- We proceed with the same criterion un- 
til we exhaust all the numbers of the artificial sequence. 
It is evident that the adoption of this procedure assigns 
to the artificial data a time order reflecting the complex 
dynamics illustrated by Figs. 2 and 3. 

At this stage, we evaluate the corresponding DE curve 
and we compare it to the DE curve generated by the 
non-shuffled real data. As earlier mentioned, the sorted 
artificial data are the same as those used to produce the 
excellent fitting of the DE curves derived from the shuf- 
fled original data. Thus, the fitting parameters are the 
same as those used for Figs. 11. We illustrate the new 
result in Figs. 12, which show that the fitting accuracy is 
as good as (and for rule No. 1 even slightly better than) 
the fitting of Figs. 11. This is a very remarkable result 
since Figs. 7 and 8 show that shuffling the data produces 
a significant effect. Thus, Figs. 11 and 12 prove that the 
memory of the data is totally under our control. 



C. Third rule in action. 

According to Lepreti, Carbone and Veltri [^3| the wait- 
ing time distribution ip( T ) is already Levy. This would 
imply that the adoption of the third rule yields an in- 
finitely fast transition from dynamics to thermodynam- 
ics. This is so because Levy distribution is stable and the 
convolution between two distinct Levy distributions is a 
Levy distribution |16|. According to our analysis, i/ , ( r ) is 
a shifted inverse power law. It is plausible that the differ- 
ence between the shifted power law distribution of Fig. 
5 and the Levy distribution of Ref. |23| is small. Con- 
sequently, the transition to thermodynamics is expected 
to be very fast. This expectation is confirmed by the nu- 
merical results illustrated in Fig. 13. The transition to 
the scaling regime is so fast that it is possible to detect a 
wide regime of linear dependence of the entropy on log (I), 
which allows us to derive for fj, the value /i = 2.138, in to- 
tal agreement with the conclusion of the earlier analysis 
done by means of rules No. 1 and No. 2. We see that in 
this case the memory of the non-shuffled data yields a S 
slightly larger than the scaling parameter of the shuffled 
data. The adoption of rule No. 3 implies a statistical 
accuracy smaller than that of the other two rules, due 
to fact there is no limitation to the jumps intensities, 
thereby decreasing the number of particles located in the 
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same cell. This has the effect of making the evaluation 
of pi and consequently that of the entropy less accurate. 
However, this disadvantage is widely compensated by the 
emergence of a much more extended scaling region that 
yields as a result a value of /i fully confirming that of the 
other two rules. 



VIII. CONCLUDING REMARKS 

We see that the uncertainty on the value of /i for solar 
flares has been significantly reduced. The current litera- 
ture, if we give the same credit to all the authors, yields 
values of \x ranging from 3 to 1.7. We provide the com- 
pelling conclusion that /i = 2.138 ±0.01. However, this is 
not the main result of his paper. We think that this pa- 
per shows that the DE method is a remarkably accurate 
technique of analysis that goes much beyond the direct 
evaluation of the waiting time distribution iJj(t). This is 
so because complex processes are characterized by two 
different kinds of memory. The memory of first kind is 
the main object of the research work done in the field of 
the Science of Complexity. To make clear the nature of 
this kind of memory, let us recall J24|] that a Markov mas- 
ter equation, namely a stochastic process without mem- 
ory, is characterized by a waiting time distribution ip(r) 
with an exponential form, thereby implying memory for a 
marked deviation from the exponential condition. This is 
why the search for an inverse power law distribution with 
a finite value of (i (the exponential distribution means 
fi = oo) can be interpreted as a search for memory. This 
is the memory of the first kind, to which the prescrip- 
tions of Ref. H are referred to. For real data, in addition 
to this form of memory, another type of memory might 
be present, denoted by us as memory of the second type, 
under the form of correlation among the values Tj. In 
this paper we have seen that this second form of memory 
is given, in this case, by the 11-year solar periodicity. It 
is possible that this form of additional memory is present 
in many other complex processes for different reasons. It 
is also evident that it is difficult, or perhaps impossible 
to reveal this form of additional memory by means of the 
direct evaluation of ip{r). This paper proves that joint 
use of the direct evaluation of ip(r) (or of ^(t)) and of 
the DE method is a very useful supplement to the ordi- 
nary technique, and that it can be profitably used to shed 
light on the dynamics behind the time series generated 
by complex processes. 

This paper yields a convincing conclusion concern- 
ing the distinction between two possible forms of non- 
stationary behavior. As pointed out in Section III, the 
claim that the waiting time distribution iJj(t) has the 
form of Eq. ( pi] ) is equivalent to assuming that the dy- 
namics of the flaring process is driven by the model of 
Eq. (12) with the assumption that the trajectories are 
injected back randomly. This is a stationary model that 
in the case where z > 2, (/x > 2), would be incompati- 



ble with the existence of an invariant distribution |?|] and 
consequently with "thermodynamic equilibrium" . The 
inaccuracy of the analyses done by the earlier work in 
this field would prevent us from distinguishing this form 
of non-stationary behavior from a genuinely form of non- 
stationary behavior. By genuinely non-stationary behav- 
ior, we mean the existence of rules changing with time. 
This form of genuinely non-stationary behavior might be 
modelled, for instance, by assuming that the parameter A 
of Eq. (TL2) is time dependent. If we make the assumption 
that the time dependence of A has a period of 11 years, 
and we make our analysis over a period of time that is 
not much larger than this time period, as we have done, 
then the process must be perceived as being genuinely 
non stationary. Our analysis is so accurate as to rule out 
the former form of non-stationary behavior and to detect 
significant effects stemming from the latter, or, equiva- 
lently, from the existence of the memory of the second 
type. 

In this paper we do not take side with cither the pro- 
ponents of self-organized criticality |2^] or with those of 
turbulence |3pl]. The dynamical model of Section III is 
inspired by the models of turbulence, but we mainly use 
it to generate artificial sequences mimicking the real ones 
with no claim that it is an exhaustive picture of the dy- 
namics behind solar flares. The fitting of Fig. 6 seems as 
good as the fitting of Fig. 1 of Ref. [^3). However, our 
analysis does not rest only on the waiting time distribu- 
tion. In a very recent paper Wheatland [ p6[ criticized 
the work of Ref. [^3| as being based on the assumption 
that rate of solar flares is constant. This is not so, as 
shown by Fig. 3. On the other hand, modelling the time 
dependence of this rate is not easy, since it does not cor- 
respond only to a 11-year periodic motion but to a much 
more complex condition, as illustrated in Fig. 4. In fact, 
this figure shows that there are many other components 
in action. This is the reason why we decided to mimic 
the time dependence of the solar flare rate sorting the 
artificial sequence in the way described in Section VII B. 
We found that this yields a fitting with the real data as 
good as the fitting between the DE curve produced by 
the artificial sequence, with no sorting induced memory, 
and the DE curve produced by the shuffled real data. 
This is, in our opinion, a strong indication that the value 
of [i = 2.138 is a genuine property of real data. On the 
other hand, the dynamical model of Section III can also 
be adapted to reproducing the modulated Poisson pro- 
cess advocated by Wheatland j[6) . This is left as a sub- 
ject of future investigation. Even in this case, the role of 
the DE analysis is expected to be crucial, and the result 
is expected to strongly depend on wether the modulation 
process involves randomness or only quasi-periodicity. 

In our notation, the power index found by the authors 
of Ref. [^U is fi = 2.38, a value that turns out to be 
compatible with the uncertainty interval associated to 
the determination of \x by means of the direct evalua- 
tion of ip{r). Our analysis establishes a connection with 
Levy statistics, in accordance again with the conclusions 
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of Ref. (23). However, we adopt a pers pec tive that is dif- 
ferent from that of the authors of Ref. 23 1. Our diffusion 
process reaches the Levy regime after the process of tran- 
sition from dynamics to thermodynamics that has been 
discussed in detail in the earlier sections. This process 
is very fast if the rule No. 3 is adopted, but it is not 
infinitely fast as in the perspective of the authors of Ref. 
|p3| who assume the waiting time distribution "0( r ) to 
obey already the Levy statistics. It is worth pointing out 
that the perspective adopted in this paper makes it possi- 
ble to take into account the time dependence of the solar 
flare rate. We do not rule out the possibility that ip(r) 
is a stretched exponential |19|. In fact, a stretched ex- 
ponential would not conflict with the attainment of Levy 
statistics in the long-time limit of the diffusion process. 
Although a truncation of ipij) at large values of r gener- 
ates a finite second moment, and consequently Gaussian 
statistics in the long-time limit, the transition to the con- 
ventional thermodynamic regime is ultra slow |p(if . It is 
known that a much earlier transition to Levy statis- 
tics occurs and that the Levy regime lasts for a very 
extended period of time. The transition to the Gaussian 
regime probably takes place at times much larger than 
the saturation time, and might be made visible only in 
the ideal case of infinitely large sequences. 
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FIG. 2. The original sequence of the solar flares waiting times. Note the logarithmic scale of ordinates. 
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FIG. 3. Number of solar flares per month from April 1991 to May 2000. 
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FIG. 5. The waiting time distribution ip(r) as a function of r. The crosses refer to real data. The dashed line is the fitting 
function of Eq. @ with A x = 31006, T = 8787 and fi = 2.12. 
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FIG. 6. f (t) as a function of r. The crosses refer to real data, and the dashed line denotes the fitting function of Eq. 
with A 2 = 30567, T = 8422 and fi = 2.144. 
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FIG. 7. DE as a function of time according to rule No. 1. The dotted straight line illustrates the slope of entropy increase 
corresponding to pi = 2.144, and S = 0.874, which is the best value of jj, afforded by the analysis of Section V. The dashed line 
is the DE curve generated by the non-shuffled real data. The solid line is the DE curve generated by the shuffled real data. 
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FIG. 8. DE as a function of time according to rule No. 2. The dotted straight line illustrates the slope of entropy increase 
corresponding to /i = 2.144, 8 — 0.5, which is the best value of /i afforded by the analysis of Section V. The dashed line is the 
DE curve generated by the non-shuffled real data. The solid line is the DE curve generated by the shuffled real data. 
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FIG. 9. DE as a function of time according to the rule No. 1. The two solid curves denote the DE curve corresponding 
to the shuffled real data, (a) The vertical bars indicate the changes of the DE curves resulting from the artificial sequences 
described in the text with T — 8422 and /x moving in the interval [2.094, 2.194]. (b) The vertical bars indicate the changes 
of the DE curves resulting from artificial sequences described in the text with n = 2.144, and T moving in the interval [7922, 
8922]. 



19 



4 




O 1 ■ ' ■ ' — ■ ' ' ' — ■ ' ' ' — ■ ■ ■ ■ 1 

1 10 100 1000 

(a) t (hours) 




1 10 100 1000 

(b) t (hours) 

FIG. 10. DE as a function of time according to the rule No. 2. The two solid curves denote the DE curve corresponding 
to the shuffled real data, (a) The vertical bars indicate the changes of the DE curves resulting from the artificial sequences 
described in the text with T — 8422 and /x moving in the interval [2.094, 2.294]. (b) The vertical bars indicate the changes 
of the DE curves resulting from artificial sequences described in the text with n = 2.144, and T moving in the interval [7922, 
9922]. 
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FIG. 11. DE as a function of time. The solid lines denote the DE curve generated by the shuffled real data, and the dashed 
lines, which almost coincide with the solid lines, denote the DE curves resulting from the artificial sequence with [i — 2.138 
and T = 8422. (a) Rule No. 1. (b) Rule No. 2. 
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FIG. 12. DE as a function of time. The solid lines denote the DE curve generated by the unshuffled real data, and the dashed 
lines, which almost coincide with the solid lines, denote the DE curves resulting from the artificial sequence with fx = 2.138 
and T = 8422 with a modulation mimicking the influence of the 11-years solar cycle, (a) Rule No. 1. (b) Rule No. 2. 
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FIG. 13. DE as a function of time, according to rule No. 3. The solid lines denote the DE curve generated by the shuffled 
real data. The dotted straight line illustrates the slope of entropy increase, S = 0.879, whith corresponds to /i = 2.138. The 
dashed line denotes the DE curve resulting from the unshuffied real data. Note the superdiffusion of the unshuffled real data 
DE due to the memory in the original signal. 
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